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We consider two weakly coupled Richardson models to study the formation of a relative phase 
and the Josephson dynamics between two mesoscopic attractively interacting fermionic systems: our 
results apply to superconducting properties of coupled ultrasmall metallic grains and to Cooper- 
pairing superfluidity in neutral systems with a finite number of fermions. We discuss how a definite 
relative phase between the two systems emerges and how it can be conveniently extracted from the 
many-body wavefunction: we find that a definite relative phase difference emerges even for very 
small numbers of pairs (~ 10). The Josephson dynamics and the current-phase characteristics are 
then investigated, showing that the critical current has a maximum at the BCS-BEC crossover. For 
the considered initial conditions a two-state model gives a good description of the dynamics and of 
the current-phase characteristics. 

I. INTRODUCTION 

A major issue in mesoscopic physics is the study of the sample sizes for which macroscopic properties emerge in 
finite systems [Tj. A typical context for such a study is provided by systems exhibiting quantum coherence, e.g. 
superconductivity or superfluidity, in samples with restricted size or number of particles: the general question is 
then to determine when and how quantum coherence takes place. The subsequent study is relevant in a number of 
situations, ranging from the investigation of superfluidity of He droplets [2] to Bose-Einstein condensation in atomic 
gases with small number of particles [3] . 

A prototypical example of these studies of quantum coherence in mesoscopic interacting systems is given by the 
investigation of the limit size of a metallic grain needed for the occurrence of superconductivity [I] : this and related 
questions are conveniently studied by using the Richardson model (RM) |51 [5J. The RM describes a system of 
attractively interacting fermions and is paradigmatic in characterizing pairing in systems with a finite number of 
fermions (5]: its relevance is also due to the remarkable feature of being exactly solvable [7i and to the fact that it 
is possible to derive the thermodynamic limit of its exact solution and show that it precisely coincides with the BCS 
solution [H1IIJ]. 

The RM has been first studied in the context of nuclear physics |7J [TU] , where the attraction leading to the pairing is 
due to the short-range nature of the effective nucleon-nucleon interaction [IF. It was subsequently shown to be deeply 
connected with the exactly solvable Gaudin magnets |12j . through the relation between the respective integrals of 
motion |13| . Using such relation, the Richardson model was then extended to more general classes of exactly solvable 
pairing-like models |14fH5] , 

The RM is particularly relevant for the study of finite-size scaling effects in the BCS theory of superconductivity 
[T§1 - |2"7] . The reason is that the classic BCS approach to superconductivity [25] in the presence of a pairing interac- 
tion violates particle number conservation [3J: number fluctuations are negligible in the thermodynamic limit, but 
important for small number of particles |5j. For this reason, the RM is used in the analysis of ultrasmall metallic and 
superconducting nanograins [5] : experiments on such systems are actually performed at a fixed number of electrons 
|29| due to their large charging energy. The Richardson model was successfully used in clarifying many features of the 
tunneling spectra of Al nanograins [29 31 , where, for instance, the spectroscopic gap between grains with an odd or 
even number of electrons was explained with the existence of pairing correlations among these |32| . 

It is a known general fact that when two superconducting or superfluid systems are weakly coupled a supercurrent 
flows between the two systems, with the current depending on the relative phase between the two superconducting 
or superfluid systems [33, J3]. The importance of this Josephson effect stems from the fact that it describes coherent 
tunneling between superfluid/superconducting systems, and this description is in most cases independent on the 
details of the microscopic description of the uncoupled systems and of the concrete physical realization of the weak 
link between them. In this paper we intend to investigate how a definite relative phase emerges between two mesoscopic 
finite-size attractively interacting systems modelled by RMs and how it is possible to extract it from the time-dependent 
many-body wavefunction: we find that this happens even for very small total number of pairs (^8— 10). This occurs 
when the "bulk" interaction (i.e., the paring interaction of the uncoupled systems) is such that the equilibrium 
properties of the uncoupled models are rather well approximated by the l&ige-N BCS theory. Once the phase is 
formed and extracted from the many-body wavefunction, it is then possible to determine the current-phase portrait 
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and study the Josephson effect in such mesoscopic weakly coupled fermionic systems. 

Our results can be primarily applied to weakly coupled ultrasmall metallic grains [JJj, but they could be also 
useful in connection with cold atom experimental setups in which the trapping potential contains a small number of 
fermions (like [35 ) and such traps are set at a distance that allows tunneling: this would be the atomic counterpart 
of superconducting ultrasmall grains coupled by tunneling terms. The fate of the Josephson effect between small 
superconducting grains was investigated to some extent in [36J , studying the dependence of the Josephson energy as 
function of the level spacing and focusing on a parameter regime where the notion of a superconducting phase variable 
is not valid. 

Another application of the Richardson model is to the study of the BCS-BEC crossover in finite size fermionic 
systems [37]. The BCS-BEC crossover is a subject which has been thoroughly investigated, also in connection to 
experimental realizations with ultracold fermions [38-40 . Increasing the (bare) attractive interaction among fermions, 
the chemical potential decreases with respect to the non-interacting Fermi energy value, so that a crossover between 
a BCS state, characterized by loosely correlated, widely overlapping Cooper pairs, to a Bose-Einstein condensate 
(BEC), in which pairs are tightly bound and minimally overlapping, can be identified |38l440] . Within the formalism 
of the Richardson model, the corresponding finite- N version of the BCS-BEC crossover [37], as well as the Josephson 
effect, can be studied. 

In this paper we numerically investigate the Josephson dynamics of two weakly coupled Richardson Hamiltonians. 
Our motivation for such an investigation is three-fold. From one side we are interested in characterizing the superfluid 
behavior of the system at finite number of particles, with regard to its phase coherence, and in investigating for which 
values of the number of particles a definite relative phase between the two systems is formed: we find that the system 
behaves coherently even for a rather small total number of pairs (as low as ~ 8 — 10) We introduce and discuss a way 
to extract from the many-body wavefunction the relative phase and its variance, so to quantify in a precise manner 
whether a well definite relative phase emerge. 

When the relative phase is well defined, we are then interested in understanding and characterizing the effects of 
the pairing interaction coefficient g (giving rise in the uncoupled models to the BCS-BEC crossover) on the coupled 
dynamics while varying the pairing interaction coefficient. We will be mostly interested to values of coupling g such 
that the uncoupled models have level occupation amplitudes close to the large- N results. 

Finally, our work aims at providing the exact Josephson dynamics between two weakly coupled Fermi systems with 
small number of fermions across the BCS-BEC crossover. Theoretical studies of tunneling of ultracold fermions across 
the BCS-BEC crossover recently appeared [411441?] . In [42] the tunneling across a barrier potential was studied by 
solving numerically the Bogoliubov-de Gennes equations at zero temperature: the Josephson current was found to 
be enhanced around the unitary limit. For vanishing barriers (i.e. large coupling between the two Fermi systems), 
the critical current approaches the Landau limiting value [42_ . Results obtained from the numerical solution of the 
Bogoliubov-de Gennes equations were compared with the analytical predictions derived from a hydrodynamic scheme, 
in the local density approximation |46| : whenever such approximation is valid (small and intermediate barriers), good 
agreement was found. In general, it is instead more difficult to obtain solutions of the Bogoliubov-de Gennes equations 
for very large barriers |47| i.e., when the coupling between the two Fermi systems is weak. Furthermore, one would 
also like to explore the exact tunneling dynamics and eventually compare it with the time-dependent solution of 
the Bogoliubov-de Gennes equations, which has been successfully used to study the dynamics of soliton solutions in 
trapped superfluid Fermi gases [50] . 

The model which is studied in the present paper, although necessarily restricted to small number of particles, 
exploits the integrability of the two uncoupled Richardson systems and allows to compute the exact dynamics when 
a state with non- vanishing initial number imbalance and/or relative phase is prepared, offering the opportunity to 
extract the dynamical phase portrait. Another advantage is that it is possible to investigate, in a simplified setting, 
how the Josephson energy depends on the interaction and the tunneling strength: we find that the Josephson energy 
has a maximum around the unitary limit, in agreement with results in literature obtained at T — in the large- N 
limit for small and intermediate barriers [421 144] . 

The plan of the paper is the following: in Section [IT] we review the main properties of a single (i.e., uncoupled) 
Richardson model. The model with two coupled Richardson Hamiltonians is introduced in Secti on |III[ where we also 
discuss the main properties of the spectrum. The Josephson dynamics is studied in Sections |IV]|V| in Section |IV| 
we introduce the considered initial states for the dynamics and we discuss the emergence of a definite relative phase 
among the two Richardson systems. In Section [V] we discuss the dynamical phase portrait, plotting the trajectories 
in the space of the relative phase and the population imbalance, and we present our results for the critical current as 
a function of the coupling. We draw our conclusions in Section |VI| 
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II. THE RICHARDSON MODEL 



The Richardson Hamiltonian is written in terms of the operators c aa destroying fermionic particles in the energy 
levels a — 1, . . . , N with spin a =f , J,: 



N N 



H = E e « (4t c «t + c U c «i) - 2 # E c U4,i c Pi c Pt- (!) 

a=l a, fi=l 

In Eq. ([!]) the the single-particles energies of the N levels, g is an interaction coefficient with the dimensions of 

an energy: in the following g is assumed to be positive (corresponding to attraction among fermions) and it models 
the matrix element of the scattering among Cooper pairs of spin-reversed fermions. The model is integrable for any 
choice of the set of energies e a - in the following we will consider them to equally spaced, according 

e a = ad, (2) 

where a = 1, . . . , N and d is the level spacing: this is indeed the choice usually done in order to recover the BCS 
physics in the thermodynamic limit (see more details below) 0[5]. 

The Richardson Hamiltonian conserves the number of fermions and, separately, of fermion pairs (doubly occupied 
levels). An essential feature of the spectrum is the so-called "blocking" effect |Sj: the states which are singly occupied, 
i.e. those in which there is only one electron with either f or 4. spin, are unaffected by the interaction and the net 
effect arising from their presence is that of "blocking" the level by preventing the scattering of the other pairs on it. 
The full Hilbcrt space is then divided into sectors with a given number of unpaired fermions and in each of these 
subspaces the Hamiltonian ([ll only couples the doubly occupied ("unblocked") levels among them, while leaving 
singly-occupied levels effectively decoupled from the dynamics. Denoting the number of pairs by M, it is customary 
to write a reduced Hamiltonian for the Nf — 2 AT paired fermions in the N unblocked levels as: 

JV N 

H = 2j2^ablb a -2g £ b%, (3) 

ct = l a, /3=1 

where we introduced the (hardcore) pair creation and annihilation operators 

b l = 4xt c U ' ba = c »ic at . (4) 

Notice that the Hilbert spaces on which the Hamiltonian ^ acts are subspaces of the full space of , characterized 
by given configurations of blocked levels. 

An Hamiltonian equivalent to ^ can be written by introducing the Anderson pseudospin-1/2 operators |51j : 
S~ = b a , = fejj, 2S^ = c^Ca^- + c^Cai — 1. In terms of the su(2) algebra generators, up to a constant, one has 

N N 

H = 2j2z«S*-2g S+Sp. (5) 

a = l a,fi=l 

Explicit solutions of the dynamics generated by the Hamiltonians ^ and §5§ has been presented and discussed in 

The exact (not normalized) eigenstates of |3| are constructed by applying a set of generalized creation operators 
B on the reference state |0) as follows: 

M 

iw)=nsK)io), (6) 

3=1 

where the reference state is the one in which no hardcore bosons are present: 

6 a |0)=0 (a = l,2,...,JV). (7) 
The explicit form of the creation operators is 

N bi 



w 

a—1 
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and the set of complex number {wj} (referred to as rapidities) satisfies the set of algebraic equations 



N 



N 



n — ' 7/1- — p_ — ' -■• 



y ct=l 3 



(9) 



The number of rapidities corresponds to the number of Cooper pairs in the state and the action of the operator ([8| 
is that of creating a boson, with a given amplitude on each level a, as results from the interaction with all the other 
Cooper pairs, which in turn is encoded in the system ([£]). In the limit g — > all the roots of the Richardson equations 
^ are real and coincide with a given subset of fields, so that each boson is localized on a definite energy level. On 
the other hand, when g is moved to nonzero values, roots can be present in complex conjugated pairs. In particular, 
for g — > 0, the ground state for a given number M of pairs is the one in which the lowest M levels are filled; in the 
strong coupling limit g — ¥ oo all the roots of this state come in complex pairs (except for the most negative one, 
when M is odd) and their absolute value diverges. The BCS equations can be obtained from this solution in the limit 
TV — > oo while keeping constant filling M/N, energy range Nd and effective coupling strength gN. In this limit, the 
root configuration associated to the ground state assumes the shape of an arch in the complex plane, whose extrema 
are at 



H±i Abcs 



(10) 



As shown in [5J EJ [HJ the link of the finite-iV results with the large- N BCS theory is provided by the fact that 
Abcs an d (J- satisfy in the scaling limit previously defined the BCS equations 



2M = Wl 

a \ 



V(e Q - M) 2 + A 



2 

BCS , 



and 



? — fj) 2 + A 



(11) 



(12) 



BCS 



The Richardson mode is integrable by means of algebraic Bethe ansatz |13l [Tf] : not only the spectrum and the 
eigenstates, but also matrix elements [53] and correlation functions [T5I [541 [55] are exactly computable. In 
particular, given two states \{v}) and \{w}) defined as in |6| with M rapidities, one can make use of 



M 



m\blb a --\{v}) = -n 



Wl 



^ vi - e a n£tj(uk - Vj)(wj - w k ) 



dct 



H - 2P a 



(13) 



where H is a M X M matrix defined as 



U?ii( w i ~ v k) 1 1 

[wj -v k ) 2 \g v. 



- 1 

Z — t o|, _ p _ Z — t 



Vk - wi 



(14) 



P is given by 



Moreover, the following relation will be also used: 
{{v}\b a \{w}) = {{w}\bi\{v}) 



j,k (Wj - h a ) 

Hili ( w i - k) 



ActH- 



& - h <*) n^fc k - v j) ( w j - w k) ' 

in which the state \{v}} has now M — 1 rapidities and the MxM matrix H~ is defined as: 



H 



3,k 



(Wj-Vk) 2 

1 

(wj-h a y 2 



-T N 1 + V 



, ■ ) k < M 

k = M, 



(15) 



(16) 



(17) 
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Figure 1: Chemical potential ll per level versus g/N, as computed from Eqs. (11|-(12| for M = N/2 — 1. Here and in the 



captions of the following figures the pairing coefficient g and the energies are expressed in units of d. 



A. BCS-BEC crossover in the Richardson model 



The Richardson model exhibits two types of crossover behavior: first, the crossover from bulk to few fermions, i.e. 
from large to small N [5]. In this case the Richardson model is used to study the corrections to the large- N BCS 
theory [9 and in general how the physical quantities are modified when the number N is not large and the energy 
scale d explicitly plays a role. Since we will numerically study the spectrum and the dynamics of coupled Richardson 
models, the size of the considered systems will be necessarily finite. Moreover, we will need to solve the equations ^ 
to determine the eigenstates, which is best done when the spacing d of the levels is kept finite while increasing the 
number of levels. It is then convenient to define an intensive Richardson gap [25 , which is related to the BCS gap 
Abcs by 

Abcs = NA (18) 

in which the quantity Abcs can be extracted from the ground state solution of the Richardson equations ([£]): it is 
found that Abcs ~ Ng [S]. In the large- N limit, the correlation functions are given by 

( b l b a) = V 2 a , (Ma) = u l , (babp) = UaVaUpVp (a ^ (3) (19) 

where the u,v's enter the BCS variational ansatz for the ground-state \GS) — Y[ a ( u a + v oMV) |0) (3J- The study of 



the comparison between the correlation functions given by ( 19 1 with those directly from the Richardson model shows 
that increasing g the agreement becomes better and better: e.g., as one can sees from Fig. 5 of [25] one has a rather 
good agreement already for N ~ 10 for g > 0.3. We can then conclude that for values of N considered in the rest of 
the paper one has for uncoupled systems a rather good agreement with large- N results. 



The behavior of ll, the real value of the extremes (10 1 of the arch formed by the Bethe roots in the complex plane 
for large values of N and g, depends in general on the filling, and it is fi cx —g for fixed values of the initial population 
imbalance, below half filling. An important point to be stressed is that in the thermodynamic limit the quantity li, 
as defined from the root configuration, tends to the chemical potential obtained for attractively interacting fermions 
in the BCS-BEC crossover [37] . 

We then can argue that the other crossover taking place in the Richardson model is the BCS-BEC one: for large N 



the parameters Abcs and fi satisfy Eqs. ( 11 1-( 12 ). Since the chemical potential changes sign for g larger than a critical 
value, therefore a BCS-BEC crossover takes place [40] . A description of the BCS-BEC crossover in the framework of 
the integrable Richardson model was given in [37], where the model (J3j was considered in the thermodynamic limit 
and it was argued there that root configurations at strong enough coupling can be used to identify the boundaries 
of the crossover. In Figure [I] we plot \i as a function of g for different values of N as determined from Eq. (Ill: for 



the considered values of N one sees that fi changes sign for g/N ~ 0.25d for M close to N/2 (notice that exactly at 
half-filling \i does not change sign). Note that, whenever M < N/2, the chemical potential becomes more and more 
negative while increasing g: at some point, it crosses the real axis to negative values, signaling the crossover. Notice 
that in the BCS scaling [5_, in which the level spacing goes to zero as the inverse of the size, the crossing point tends 
to a finite value of g in the thermodynamic limit, whereas in the considered equally spaced model Q the crossing 
occurs at a value of g which is instead cx N . 
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III. COUPLED RICHARDSON HAMILTONIANS 



In this Section we introduce the model studied in the rest of the paper featuring two Richardson models coupled 
by a tunneling term |33[ [34] : 

H = H R + H L + H t , (20) 

where H R and Hl are the "right" and "left" Richardson Hamiltonian, written in terms of the right and left operators 
b a ,R.,b a ,L (the fermionic operators will be denoted by c aaiR and c„ ffi i with a =f,J,). We consider the simpler setting 
in which the two models have the same value of the coupling g and the ss,me energy levels £ a \ 

N N 

H S = 2J2 e«bl s b a , s -2gJ2 b lM* ( s = L > R l ( 21 ) 

Q = l a, 0=1 

with the e a 's equally spaced and given by ^ . The number of levels is taken to be equal to N both for the left and 
the right systems. The total number of pairs in the system is denoted by Mt - we will also denote by Ml and Mr 
the operators of the number of pairs in the left and right system: M s — J2a=i b a s^a.s (with s = L,R). 

We write the tunneling term describing the hopping of a single fermion from one system to the other in the form 

N 

H t = -vJ2 E { c L,l^,r + h-c) (22) 

u=\,i a, ,3=1 

(with rj having the dimension of an energy). Following the usual approach initially introduced by Josephson [33 , 
using second-order perturbation theory one can derive an effective Hamiltonian for small values of A (corresponding 
to the regime of weakly coupled Richardson Hamiltonians): it turns out the this effective Hamiltonian can be written 
only in terms of the pair operators [36J, greatly simplifying the study of the dynamics. 

Since the uncoupled Hamiltonians (21 1 contain only interactions among pairs, the eigenstates of ([!]) can be classified 



in terms of their seniority v, i.e.. the number of the unpaired electrons. The second order effective tunneling term 
can be written as: 

H(2) =-E E E gT |QL/W) {ahPRa ' M H T , (23) 

in which the sum runs over all the possible intermediate states that can be reached from a z/-seniority couple of states 
\N/2 + v) L ® \N/2 + v) R , by removing an electron of spin a from the level /3r of the right grain and adding it on the 



level «i on the left grain (or viceversa). In ( 23 1 the quantity E aL p RV is the corresponding excitation energy relative 
to the initial state. 

Following [36J, it is possible to limit the space of states on which the intermediate sum runs over to the lowest 



energy ones, when acting with (23) on the lowest-energy states of the two uncoupled systems in which all fermions 
are bound into Cooper pairs. In facts, the energy E aL p RV includes the energy needed to break a pair and the effect 
of the blocking of the states on the collective excitations on both subsystems. 

Across the whole BCS-BEC crossover, the breaking of a pair associated with the tunneling of a single electron is 
energetically depressed: processes like the ones depicted in Figure [2^ a)-(b) are suppressed by a factor 1/ Abcs m the 
dynamics, since they involve both the breaking of a pair, with an energy cost equal to the BCS gap Abcs and the 
blocking of a level, which affects all the levels and has therefore an energy cost roughly proportional to N. At second 
order in the fermion tunneling, it is more convenient to reach a final state in which only Cooper pairs are present 
|33| . Moreover, single-fermion tunneling does not produce a current in the absence of an applied driving force, so 
they will not affect the current. This is true in particular for processes like the one in Figure [2jc), which reproduces 
the initial state and can be included in a redefinition of the energies of the unperturbed system. We are therefore led 
to consider as dominant the coherent pair tunneling, which involves both the electrons of a Cooper pair and can be 
directly written in terms of the bosonic operators b* a L , bp t ji or 6 q .l6^ r , as in Figure J^d) . Assuming the two systems 



to have a definite relative phase (as it will be checked and discussed in Section IV I , the coherent tunneling involves 
a phase shift on the state in which it takes place and a corresponding variation of the relative number of fermions 
SN f = ±2 [see Figure^ 



We therefore focus on coherent pair tunneling, for which the effective Hamiltonian is 



* ti a + hp 

a,p H 
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Figure 2: Second-order processes associated to the Hamiltonian |23 



where E a = + A^ cs , £ Q = e Q — /i and A = 2r; 2 / Abcs- As can be seen in Eq. (24 1, we decided to scale the 
tunneling coefficient with Abcs since this is the relevant scale throughout the crossover and it is a finite quantity in 
the thermodynamic limit; furthermore, the form (24) ensures that the tunneling acts as a perturbation also on the 
BCS side and for small values of g. 

The form d24J) is particularly relevant since it formalizes the fact that preparing the system in its ground state and 



adding a weak fcrmionic tunneling term to the uncoupled Richardson Hamiltonians does not destroy the Cooper pairs 
picture: this provides a major simplification in the problem, allowing for to study the Josephson problem only in 
terms of hardcore bosons since the subspaces with different seniority will not be accessed neither by the "single-site" 
dynamics of the uncoupled Richardson systems, nor by the coupling between different sites (see a discussion on the 
single- fermion tunneling effects at the end of Section IV ) . 

The ground-state state of Hamiltonian H = Hl + Hr + was studied in and the behavior of the Josephson 
energy investigated as a function of the level spacing. Integrable versions of coupled pairing Hamiltonians was 
proposed and studied in [56 - 59J , while an analysis of the spectrum of two weakly coupled Richardson Hamiltonians 

with H T cc J2 a .p (f a ,L b P,R + & ",£ & /3,r) was presented in [Ml [50]. 

In the following we consider the Hamiltonian H = Hl + Hr + H^ 2 \ with given by (24 1 and we investigate 
the properties of its spectrum and the dynamics starting from an state at time t = having an initial relative phase 
difference and/or an initial population imbalance: we are interested to ascertain for what values of N a relative phase 
difference 5(f> is well defined, and to study the dynamics in terms of the time evolution of S<f>(t) and SM(t), where SM 
is the difference between the number of pairs of the two systems defined by Eq.(31l. Note that, in general, the gap 
and the chemical potential appearing in (24 1 will be functions of time as well. However, for the sake of simplicity, 
we will consider them as constant, which in the present case stands as an approximation valid for small fractional 
population imbalance. 

The initial state is prepared in the following way (see more details in Section IV I: the uncoupled system (A = 0) is 
initially in the ground state, characterized by a definite occupation number on the left and on the right parts - then, 
at time t — 0, the coupling A -C d is turned on and the quantum dynamics is studied. 

Integrability plays an important role both in the study of spectrum and dynamics: it gives the exact eigenstates of 
the two uncoupled systems and, most importantly, the exact hopping matrix elements. It also provides an efficient 
truncation mechanism to select the most important eigenstates in the dynamics: as we discuss in the following, one 
can avoid to diagonalize the Hamiltonian written on a basis of the full Hilbert space of the coupled problem and 
instead limit its size by restricting only to a subset of states. 

Operatively, we start from the basis of the exact eigenstates of the two uncoupled models, with A = 0, with the 
number of pairs on the left (Ml) and on the right (Ml) separately conserved. Once fixed the total number Mr of 
pairs, the factorized basis Sm t is split into subsectors, each of them characterized by the occupation number of the 
left model Ml and that of the right one Mr, such that Ml + Mr = Mr. Denoting by Sm a basis of eigenstates of 
^ for the subspace with given number M of pairs, the fixed- number subspaces are spanned by 



>M l ,Mr 



= $ 



(Mx, 



4> 



(Mr) 



(Mr) 



€ Sm R } ) 



so that the factorized basis is 



min(Af,M T ) 

Sm t = U Sm,m t 

M=max(0,2N-M T ) 



-M ■ 



(25) 



(26) 



It is possible to show that many states in Sm are effectively not involved in the dynamics and consequently reduce 
the space of quantum states to a computationally manageable size: to see this, let first consider the limits g —¥ and 
g — > oo. In the non-interacting case, the single-level occupation numbers are good quantum numbers for the system: 
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Figure 3: Particle-hole states. 



it follows that all the excitations above the Fermi sea ground state induced by the coupling, in the regime in which 
the tunneling coupling is small (X/d <C 1), are the states in which one particle is missing from the Fermi sea or one 
particle is added above it. These are a subset of the "particle-hole" states, obtained from exciting one pair from below 
to above the Fermi level, which are instead there at second order. 

In the opposite limit g — » oo, it is useful to rewrite ^ in terms of spins, obtaining the spin Hamiltonian in the 
strong coupling limit g — > oo, the Hamiltonian ^ reads [23 



H 



-2g( 



S-S- (S z 



S 1 



(27) 



(where S — J2 a &a) an d it conserves the total spin of the state and its z-projection. Numerical solutions of the 
Richardson equations show that the rapidities can either diverge proportionally to g or remain finite, with real part 
which lies "trapped" between two energy levels. In the strong coupling limit, the tunneling Hamiltonian. Conse- 
quently, the states group into highly degenerate total spin subspaces |23| . In the strong coupling limit, the tunneling 
Hamiltonian (23) simplifies as well: the BCS gap diverges linearly with g and all the pairs of levels in (24 1 factorize 



a common term, yielding 



H (2) 



AA 



BCS 



2 

BCS 



^tot,L^tot,R 



h.c. 



(28) 



(where Sl — Y] a S a _L and where Sr = S at R). The ground state is the unique state in which all the rapidities 
diverge in the strong coupling limit and it is the one with highest (total) spin. The relation between the number r of 
diverging roots at strong coupling and the eigenvalues of the spin Hamiltonian (27 1 is r — s — s z [23 , being s(s + 1) 
and s z the total spin projection along the z axis. One then sees that it is sufficient to restrict the single-site Hilbert 
space to the root configurations with one less (or one more) rapidity and the same number of rapidities which diverges 
at large g, i.e., again the ground state of the new sector: therefore, no new state is needed. Although the previous 
arguments are valid in the two limiting regimes g — > and g — > oo, we numerically compared the results with exact 
diagonalization (for N — 6) or the effect of adding more total spin subspaces to the dynamics (for N = 8). In all the 
tests we performed, results in excellent agreement were found. 

Algorithms for connecting the number of roots that eventually diverge to the initial state configurations have been 
given in |25M61j : the included states are exemplified in Figure[3]and consist of the evolution in g of all the configurations 
in which, in the weak coupling limit, one particle is excited from the Fermi sea to right above its surface or from the 
Fermi energy to one more energetic state. 

The algorithm used for solving the Richardson equations numerically is based on the one described in [62_ . To obtain 
eigenstates at a given g, one starts from g — 0, where the rapidities that solve the Richardson equations are known 
within good approximation. It is therefore possible to solve numerically (9| for some values of g around zero: the 
coefficients of the polynomial having these rapidities as roots are computed. One then extrapolates these coefficients 
to a new value of the pairing, in steps Sg ~ O.Old, and compute the roots of the extrapolated polynomial, using them 
as a starting guess for the numerical solution of the Richardson equations. The procedure is iterated up to the desired 
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value of g [62 (see more details in [63 ). This algorithm allows to solve every configuration for sizes N < 10, which 
we used in this paper. A numerical procedure for dealing with general Gaudin models has been presented in |64| . 
Once that the factorized basis has been determined, we write the tunneling term by using Eq. (16 1, and diagonalizc 
the resulting Hamiltonian. 



A. Properties of the spectrum 

In Figure [4] we plot the energy spectrum for two coupled Richardson models as a function of the tunneling parameter 
A for three different values of g. It is seen that the effect of a weak tunneling on the spectrum depends essentially 
on the coupling among fermions: one can clearly identify a regime of nearly non-interacting particles, in which the 
quasi-degeneracy of the levels is given by the number of ways of promoting one or more particles in an excited level 
to obtain a given energy (degeneracy is a consequence of the choice of equally-spaced levels). In this regime, the 
perturbation splits the levels of one band as far as the band spacing, hence giving rise to a spectrum in which the 
original degeneracies are not seen any more. 

On the other hand, in the strong coupling regime states group into eigenstates of the total angular momentum, as 



seen from the spin representation (27 1. Since the distance among the energies of these subspaces is of order g, in this 
regime, even a tunneling term of several times the gap cannot mix the different subspaces among them. 

In the crossover region, the strong coupling subspaces are already quite defined, but not far one from the other. It 
follows that a sufficiently strong perturbation can still hybridize them. To better illustrate this point, we may evaluate 
how much the levels are shifted by turning on A: however, the absolute value of the shift should be compared with 
the level spacing in a situation where levels are well-distinguishable (intermediate couplings) and the band spacing in 
the presence of strong degeneration (g — > or g — > oo). 

A convenient way to proceed is to divide the all spectrum in a certain number of interval (let A^ in this number) 
and count how many levels lie in each interval: we define the quantity 

Xm = levels in the m th interval)(X) — levels in the m th interval)(X = 0). (29) 

The average of this quantity with respect to the interval index m is obviously zero. The relevant quantity is instead its 
standard deviation a x : the result is plotted in Figure [5] left part. The quantity <r x has a maximum around g/d ~ 1, 
corresponding to the point in which the degeneracies of the noninteracting picture are already destroyed, while the 
energy bands of the strong coupling regime are not evident yet. This result is of course related to the occurrence of a 
crossover between weakly attracting fermions and tightly bound pairs: for the uncoupled systems, from the point of 
view of the energy spectrum the crossover reflects itself on the creation of energy bands out of the pair levels, which 
are more and more separated by increasing g. This is also seen at the level of the coupled spectrum, The doublet 
structure characterizing coupled noninteracting systems is melted into an highly-degenerate band structure. 

In the right part of Figure [5] we plot the energy difference between the first excited state and the ground state in 
a system with an odd number of particles as a function of g. It turns out that, as long as the gap opens more and 
more, the energy difference between the components of the level doublet reaches a maximum splitting. 
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Figure 5: Left: standard deviation a x vs. g associated with the level structure, given by (Nun = 100, A = 0.05, N — 8, 
Mt = 6). Right: energy difference between the first excited state and the ground state as a function of g for the same values 
of the parameters. 



IV. EMERGENCE OF A DEFINITE RELATIVE PHASE 



In this Section we explain what initial states have been considered and we discuss how a definite relative phase 
emerges for small number of particles and our algorithm for obtaining it from numerical data. 

The initial state \*f?(t — 0)) is prepared in a linear combination of ground states of the uncoupled systems having 
a different number of pairs and therefore a population imbalance: this state is evolved in time with the dynamics 
generated by (20 1, with A ^ in the tunneling term. More precisely, at t = the system is generically in a linear 
superposition of two states with a given number of pairs (we choose M in a system and M — D in the other): 



the total number of pairs is conserved during the time evolution and is Mt — 2Mo — D. Denoting by $ 
lowest-energy state with M pairs of either the left or the right system, we prepare the system in the state 



(L,R) 
M 



l*(* = o)> = 



1 



4> 



(R) 

M -D 



4> 



(L) 

Mn-D 



M I 



the 



(30) 



given the limitation on the total number of pairs (Mt < 10), we will most consider D = 1, therefore creating as initial 
state by a linear superposition of Mq in a system and Mo — 1 in the other, and D — 2 (typically we choose Mo = N/2 
or M = N/2 ±1). 

Given |\&(t)), one can compute the pair population imbalance as 



SM(t) = (tf (i) \M L - M R \ *(<)> = ( (t) 



E 



*(*) 



(31) 



(the population imbalance SN f is just 6N f = 2SM). Notice that with £ = it is (*(i = 0) \M L \ $f(t = 0)) = M . 
According to the notation of |65| . we will denote by z(t) the fractional population imbalance: 



z(t) 



SM(t) 
M T 



(32) 



Using the state ( 30 1 one simply obtains 



6M(t = Q) = D 1 YT ^: 

therefore varying the parameter £ one can choose different initial population imbalances (with \5M(t = 0)| < D). 
The dynamics is then studied turning on a small perturbation (X/d — 0.01 — 0.1 in our runs) and compute the time 
evolution of the state after exact diagonalization the Hamiltonian. The main limitation of this protocol arises from 



the consume of RAM by diagonalization subroutines: by limiting subspaces appropriately, as discussed in Section III 
one can study systems up to TV = 10 levels (both on left and right systems). 

An important issue we want to address in this Section, arising from the fact that we can treat the exact quantum 
dynamics of the coupled model only for a limited number of pairs, is whether a definite relative phase emerges at small 
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sizes. In the presence of the tunneling term (23), eigenstates will in principle be written as a combination of many 
of the factorized states of the two uncoupled Hamiltonian. Nevertheless, as we will see, when the initial population 
imbalance is small, the number of involved states is rather small. Moreover, even for higher particle imbalance, 
when the tunneling is weak and the pairing strength is strong enough, the Hilbert space of each system organizes in 
subspaces, labeled by eigenvalues of the total spin (see Section [IT]) . It follows that in most cases, even if the exact 
states involved are many, the corresponding energy eigenvalues are not very different, therefore the time evolution 
takes place with nearly definite phase. 

Note that in our canonical setting the expectation value (W(£)| b a l/r \^(t)) is always vanishing, since the b's 
operators does not conserve the number of particles. However, we can define time-dependent phase differences between 
one level in a system and a level in the other system systems by the use of the formalism of Section [TT] evaluating 
the dynamical two-point functions (for the uncoupled systems, dynamical two-point correlations have been studied 

MM)- 

From the correlation function ('l'(i)l b* a ibp,R one can extract how much the phases of two distinct levels differ 

at a given time. In particular, we considered two different procedures for the choice of the levels, which can be tested 
one against the other, and define: 

w a (t) = (*(t)\bl L b a , R \*(t)) (33) 

and 

z a (t) = <*(t)| bl L b N/2tR |*(t)>. (34) 



In (34 1 the subscript refers to the level on the left system and a reference state is taken on the right system (arbitrarily 



chosen to be the level N/2); conversely, in (33 1, the level is chosen to be the same on both systems. We define a 
relative phase between levels as 

WaW^kWIe^" 1 (35) 

and 

z a (t) = \z a (t)\e iS ^ a \ (36) 

The functions 5<j) w {t; a) and S<j) z (t; a) are functions of both time and level index. It is then necessary to verify whether 
the levels have small phase difference: to do this, we define the level average 



1 N 



and their standard deviation <J w , z (t) [with <r^ z (t) = (l/N) X^=i(<5</v z(*i a ) ~ 5(j) WtZ (t)) 2 ]. The time evolution of the 
mean values 54> w ,z(t) is reported in Figure pj one sees that already for Mj- = 8, one has relatively small values of 
g where the two definitions of the relative phase are in good agreement for most of the times. The two definitions 
S(p w (t) and S<fi z (t) are expected to agree only when the two systems show coherent behavior, and the phase difference 
between them is, within a good approximation, given by the phase difference between any two levels chosen. We 



checked that choices other than ( 33 (-(34 1 give practically the same results when a relative phase is well defined. 

In order to have a definite relative phase one has to check that the average values S<j> w ^ z (t) should be (possibly for 
most of the considered times) much larger than their standard deviations a w , z (t): as shown in Figures [7] and [8] (d one 
respectively for g = 0.2d and g = 0.6d) this condition is rather well verified also for a number of pairs Mt = 8. One 
also sees that for g = 0.2d the agreement is less good, as expected also from the fact that - as discussed in Section [IT|- 
the uncoupled systems have significant deviations from the large- TV limit. We also observed for the considered values 
of g a significant degradation of the relative phase for even smaller total number of pairs, e.g. as low as Mt = 4. 

Information about the phase difference averages and their standard deviations at every given time is useful, but 
we can complement it with their averages in time: to this purpose, we consider the mean of the standard deviation 
presented above over sufficiently long times (several periods) 

c£?, = -!- / " ma ° ' <T W<z {t')dt' '. 

To establish a comparison, we need to evaluate also the mean phase difference among the condensates. This is an 
oscillating quantity, having vanishing average on time: we then compute the average of its square: 

SL = \ v^- f* m " Hi w (t')dt'- 
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Figure 6: Phase differences 5cj> w (t) (solid blue line) and 8(j> z (t) (dotted red line) vs. t for the coupled systems with N — 8 levels 
each, total number of pairs Mt = 8, pairing strength g = 0.6, tunneling parameter A = 0.1, D — 2 (corresponding to Mo = 5) 
and initial imbalance z(t = 0) = 0.25. Time here and in the following figures is in units of h/d. 




Figure 7: Phase difference means 8<f) w ^(t) (dashed blue lines) and standard deviations a w , z (red solid lines), as determined 
from the correlation functions to (left) and z (right), for two coupled grains with the same parameters as in the previous Figure. 




Figure 8: Same quantities as in Figure [TJwith g — 0.2 (other parameters unchanged). 
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Figure 9: C S J (lower dashed line), S S J (upper dashed line) and Cf^ (lower solid line), S s z ^ (upper solid line) vs. g. Parameters 
are as in Figures [6p| N = 8, M T = 8, A = 0.1, D = 2 (moreover, t max = 1000). 



In Figure 9 we plot C 5 ^ and S 5 ^ with both the definitions p3|-d34l. One sees already for g > 0.3d a very good 
agreement C^f and Cf^, and both significantly larger than the time averages Sf^ z . For this reason we are going to 
denote as 5tfi the relative phase difference, omitting the indexes w, z. One also sees that for small g the relative phase 
is not defined, as expected, since the relative phase is comparable with its variance. 

As a function of the pairing parameter g, from Figure [9] one sees that the higher is the value of g, the more the 
system shows a definite relative phase. We also observed that the smaller is the tunneling parameter, the sooner (in g) 
a definite phase is established. Similarly, a small initial imbalance allows for a definite phase to emerge for relatively 
small values of g. while - for the considered values of N - stronger pairing is necessary if states with larger initial 
imbalances are selected. This is due to the fact that the initial state is projected on few states in the lowest part of 
the spectrum when the initial population difference is small. Conversely, larger population imbalances at t = are 
projected to many states in the middle of the spectrum, each having its own energy. 



We pause here to comment about fermion tunneling: as a matter of fact, the original tunneling Hamiltonian (22 1 is 



written in terms of fermionic operators, while the result that the phase coherent behavior is established with relatively 



small pairing and/or total number of pairs is obtained with the bosonic approximation (24 1, acting on the restricted 
subspace of blocked levels. Since at small g pair-breaking excitations may play an important role, a natural question 
to ask is whether the presence of fermionic degrees of freedom, aside of bosonic pairs, may spoil the phase-coherent 
behavior of the systems for sufficiently large pairing. The issue can be rephrased into the question of whether the 
initial state, during the evolution generated by the coupled Hamiltonian, containing a fermionic tunneling term, may 
give rise to a huge number of states in which two or more electrons are not paired, evolving incoherently with respect 
to the states in which only pairs appear. 

These states have to be written as linear combinations of the factorized states of the two uncoupled Hamiltonians. 
On each site, the energy of such states can be exactly computed for any value of g. In order to have an estimation 
of a lowest bound for the energy, we can consider a state in which the most energetic pair is broken and one electron 
is promoted into the next level, which reduces the number of pairs by 1 and the number of unblocked levels by 2, as 
seen m Section [IT] The energy of the lowest pair-breaking excitation has been considered in [23] and it reads: 

E pair ~ £M + o £M+1 - g(M - 1)((N - 2) - (M - 1) + 1) (38) 



The bare energies in the first term of (38) do not depend on g, unlike the ground state energy, all the pair-conserving 
excitations and the second term in the previous equation. It follows that, by taking the pairing strength sufficiently 
high, all pair-breaking excitations can be made to lay at arbitrary energy above the ground state and are therefore 
suppressed with respect to pair-conserving excitations. 

Checking explicitly that the insertion of states with unpaired electrons does not spoil the phase relation requires 
much larger computational effort, in that the Hilbert space should be enlarged to the ( ) { m™) configurations in 
which the m electrons can "block" part of the TV levels, with fixed number M of pairs. We can therefore qualitatively 
rely on the standard argument based on the presence of a gap preventing single-fermion tunneling: note that this 
should already hold for values of g > 0.25, as previously discussed. 

We also mention that, even if the phase is quite well defined, residual fluctuations can still be observed, in such a 
way that the widest, slowest oscillations are superimposed with faster and narrower ones. We find convenient to isolate 
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the former ones by computing time averages on intervals much smaller than the period of the largest oscillations: this 
allows to better understand the structure of the dynamical diagrams discussed in the next Section. An example of 
the procedure is provided in Figure |10| 



V. PHASE PORTRAIT AND CURRENT-PHASE CHARACTERISTICS 

In this Section we first draw the population-phase dynamical portrait z(t)-6<f)(t) as it has been done for bosonic 
Josephson junctions |65l I67j and we determine the current-phase characteristics, which is a typical tool used to 
characterize the behavior of a Josephson junction |341 155] , From the solution for the quantum dynamics one can 
extract the dominant period of the population oscillations and determine the Josephson frequency. We also comment 
on the determination of a two-state model giving a good description of the dynamics and of the current-phase 
characteristics for the considered initial conditions. We observe that most of our simulations are done for the initial 



state (30 1 with D = 1 build by a linear combination of a state with Mq = N/2 pairs on a system and Mq = N/2 — 1: 
this state has a maximum value for \SM(t — 0)| equal to 1. For this initial state the relative phase difference 5cf> is well 
defined for a total number of pairs > 6 and for g > 0.3 (see Figure^]), where the expectation values for the correlation 
functions are already rather similar to the large- TV BCS findings [25]: we can then explore the crossover region (which 
is around g/d ~ 0.25N). In the final part of the Section we consider D = 2 and initial imbalance SM(t = 0) = 2: 
the phase turns yet to be again rather well defined (but a larger values of g) , but we cannot practically explore larger 
initial imbalances (i.e., larger values of D) since with our maximum value of pairs ~ 10 the relative phase is well 
defined only for very large values of g (well beyond the crossover point). 



In Figure 11 we plot the number-phase portrait where we plot as a function of time both S<f>(t) and SM(t) for 



different values of SM(t — 0). It is also possible to study the diagram while varying the initial phase in the initial 



state (30 1, as shown in Figure 12 

One sees from Figures |ll|[l2 that even for a small total number of pairs (My = 7) the phase diagram in the plane 
8<f> — z shows a remarkable agreement with a "pendulum" law of motion in the small oscillations regime when the 
initial imbalance is small. Furthermore, as the initial displacement or phase difference becomes larger, significant 
corrections are seen. 

A way to understand such results is to introduce a two-state model [34, 69 : computing the overlaps of the initial 



state (30) having D = 1 and \SM(t — 0)| < 1 with the many-body eigenfunctions of the full Hamiltonian (with A 
small), one sees that the largest overlaps are with the ground and the first excited states. Given this one expects 
that the dynamics is well explained by a simple linear two-mode model involving such two states. The dynamical 
equations of the Feynman two-state model are reviewed in Appendix [A] for the linear two-state model here considered 
the phase difference 8<j) does not overcome the value n/2, i.e., if \<f>{t = 0)| < tt/2, then \<fi(t)\ < ir/2. As seen in Figures 
TTfl2 this property is clearly observed in the numerical results (we also checked it with exact diagonalization) . The 
property is typical of the linear two-mode model and it is connected with the fact that the main contributions to the 
time-dependent wavefunction arise from the first two lowest-lying states of the interacting system with equal weights. 
We now focus on the pair current between the models: we define the current / as the time derivative of the 
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Figure 11: Dynamical phase-portrait z-8<j> for different values of the parameter £ in the BCS regime, with N — 8, Mt = 7, 
D = 1, g = 0.57, A = 0.05. The chosen values of £ are £ = 0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9 corresponding respectively to 
<5M(t = 0) = 0.98, 0.92, 0.83, 0.72, 0.60, 0.47, 0.34, 0.22, 0.10. 
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Figure 12: Phase portrait for different initial phases in the BEC regime (TV = 8, M T = 7, D = 1, SM(t = 0) = 1, g = 9.7, 
A = 0.05). 



occupation number of the left subsystem 

d 



I(t) = ^(*(t)\M L \9(t)) = -(9(t) 



*(*) 



From Eq. ( 39 1 one finds 



M(t) = i [H, M L ] = i 



H (2 \Mr 



(39) 



(40) 



As discussed in the Appendix [A] for the linear two-state model the current is proportional to the tangent of the 
phase difference: / oc tan (50 [see Eq. ( |All )]. The current-phase characteristic can be therefore written as 



I(8<t>) — I c {g, A) tani 



(41) 



and the critical current I c can be fitted from numerical data. An example is given in the left part of Figure |13[ We 
find that the critical current has a maximum around a finite value of g, as shown in the right part of Figure [13] for 
the considered values of N the maximum is at g ~ 1, close to the unitary regime. I c can be fitted in the form 



Ic(5,A) = /oA- 



-c/s 2 



g 



lo depends mostly on N. Notice that the relation (42 1 has a maximum at g* = \/2c. For the parameters of Figure 
we find c ~ 0.27, nearly independent on A. 
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Figure 14: Number-phase diagram for JV = 8, D = 2, 5M(t = 0) = 2, M T = 6, A = 0.05 and g = 0.4. 



We stress that the fit needed to identify the critical current is done using the linear two-mode model: the validity 
of the fit relies on the fact the two lowest levels are the ones mainly involved in the dynamics, which is the case for 
small imbalances (D = 1). Deviations are observed for larger values of D, as we are going to discuss. 

It is an interesting issue to explore what happens when more levels, inserted in a band structure as the one described 
in Section III A| participate the dynamics: with D — 2 and SM(t = 0) = 2 the phase diagram shows a typical ellipsoid 
form. An example of number-phase portrait is given in Figure [l4j We see that the phase range depends only on 
the interaction, while the amplitude of the population oscillations depends on the initial relative phase given to the 



system through ( 30 ) 



The numerical study of the current phase characteristics reveals that for D — 2 the relation ( 41 1 does not provide 
a good way of fitting the critical current: the numerical results are plotted in the left part of Figure [15] We find that 
a good approximation of the current-phase characteristics is given by 



7(^)=/ c ( ff ,A)sin- 



(43) 



with I c given by (42 1 [70], as it can be seen in the right part of Figure 15 We observe that such a dependence for the 
current-phase characteristics was found for a weak, point-like barrier in the WKB approximation in the Bogoliubov-de 
Gennes equation [71]. Since for large N we expect a dependence cx sin S(f> [47], we attribute the result (43) to the 



small N considered: further numerical investigations with larger number of levels are needed in order to obtain the 
current-phase characteristic for intermediate and large N for the coupled Richardson models. 
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Figure 15: Left: current-phase characteristics for g — 3.8, A = 0.05, D — 2, 5M(t — 0) = 2, N = 8, Mr = 7 - blue circles are 
numerical points obtained from the quantum dynamics, and the red line is the fit according Eq. (431. Right: critical current 
vs. g - the blue circles are the numerical results, while the red line is Eq. (|42[). 



A check of Eq. (43 1 and of the data presented in Figure 15 can be obtained by doing the Fourier transform of 5M(t) 
with D — 2: as a function of g, to a very good approximation the dominant frequency of SM(t) (i.e., the Fourier 
component with the highest weight) turns out to be proportional to the critical current given by (42). 

An important prediction of the nonlinear two-state model is that there is a critical initial imbalance for which self- 
trapping occurs |65j : given the limitation on the maximum value of D, we cannot explore larger initial imbalances. 
What is observed, instead, is that the amplitude of the fastest oscillations of SM(t) is increased and that the period of 
the slowest ones is decreased more and more, as l/g. The time period of both SM(t) and 5cf>(t) become larger and the 
oscillations exhibited by 5<p{t) (as the ones seen in the left part of Figure 10| become as well larger. The scenario is 
that of a large crossover to a confined regime, in which the occupation oscillations have infinite period at g very large. 
This may be a finite- N effect, and one could expect that this eventually leads to a transition in the thermodynamic 
limit. 

The initial phase can also be varied with initial imbalance 8M(t = 0) = 2. It is interesting to note that the for 
most of the values of g, the phase runs: nevertheless, the time evolution of the mean phase locks it around some 
large-period oscillations. We conclude by observing that similar results are found decreasing the coupling A: further 
investigations to study self-trapping effects at very small values of the coupling are needed. An analysis of larger 
imbalances and larger N (eventually with very small coupling) is therefore needed to study self-trapping effects, and 
more in general non-linear effects, through the crossover. 



VI. CONCLUSIONS 



We have studied the emergence of a definite relative phase between ultrasmall metallic grains (and in general finite- 
size systems of attractively interacting fermions) modeled by weakly coupled Richardson models. We have introduced 
and discussed a way of extracting the relative phase and its variance from the many-body wavefunction, in order to 
precisely quantify whether a definite relative phase emerges. 

We have also related the coherent behavior to the spectrum of the coupled systems and suggested a criterion to 
characterize the crossover between the BCS and BEC regimes, showing that these regimes are clearly distinguishable 
by the spectrum of the coupled models. 

Moreover, we have performed a numerical analysis of the exact dynamics of the two weakly coupled Richardson 
Hamiltonians, after a weak tunneling term is turned on. We used a linear superposition of the eigenstates of the two 
uncoupled systems, with a different number of pairs (D being such difference), as initial states: these states are then 
evolved according to the full Hamiltonian including the tunneling Hamiltonian, weakly coupling the two systems. 
We found that a definite relative phase difference emerges even for a small numbers of pairs (~ 8 — 10). Therefore, 
the current-phase characteristics could be obtained for values of the bare pairing strength for which the equilibrium 
properties of the uncoupled models are well approximated by the BCS theory. We showed that, for small initial 
imbalances (D = 1), a two-state model gives a reasonably good description of the dynamics and of the current-phase 
characteristics. 
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Finally, we have presented the critical current as a function of the pairing parameter, finding that it has a maximum 
around the unitary regime, even with a number of pairs ~ 8. The phase portrait was studied for small initial imbalances 
(D; 2). 

The requirement of having a definite phase difference among the two systems with a limited total number of pairs 
(< 10) prevented us to analyze values of the initial population imbalance (D > 2): for these large initial imbalances 
the relative phase is well defined only for very strong pairing interaction, well beyond the unitary limit and deep in 
the BEC regime. Further numerical investigations are required to consider larger sizes and larger initial imbalances 
(eventually with very small tunneling couplings), which may generate a definite relative phase across the BCS-BEC 
crossover: it is expected that a proper finite-size scaling may be crucial to identity non-linear self-trapping effects. We 
moreover regard as interesting the investigation of the effects on the relative phase of single-fermion tunneling terms: 
these terms might give a contribution on the BCS side of the crossover and produce a degradation of the relative 
phase, which should eventually form for larger sizes. Similarly, it would be stimulating to compare (eventually for 
larger systems) the results obtained from exact dynamics with the ones obtained using time-dependent mean-field 
approaches. 

The rapid growth of the computational cost with the size of the systems represents a limitation on the total number 
of pairs as well: the Hilbert space could be further reduced in the strong coupling regime, yet not throughout the 
whole crossover. We conclude that it stands as an open issue, certainly deserving future work, how our findings scale 
with the size of the system. 

Our results can be applied to weakly coupled ultrasmall metallic grains and to cold atom experiments in which 
traps with few fermions are set at a distance that allows tunneling: the individuation of the relative phase between 
nearest neighboring sites makes possible in perspective to study Josephson dynamics and self-trapping systems also 
for larger imbalances, and to check the validity of two- and multi- mode ansatz. 

We finally observe that in this paper we focused our attention to weakly coupled Richardson models, discussing the 
formation of a relative phase and the Josephson dynamics for a class of considered initial conditions. The extension of 
our method of defining a relative phase to the problem of the formation of a relative phase between general interacting 
(both integrable and non-integrable) mesoscopic systems could be relevant in a rather broad class of physical systems, 
including weakly coupled ultracold finite Bose gases, and it is in our opinion an interesting problem, worthwhile of 
future studies. 



Acknowledgments 

Discussions with G. Sierra, A. De Luca, T. Macri, A. Smerzi, L. Amico, R. Scott, L. Pitaevskii and S. Stringari are 
very gratefully acknowledged. F.B. also thanks A. De Luca for collaboration on the implementation of the numerical 
solution of the Richardson equations. 



Appendix A: Dynamical equations for the two-state model 

A general description of the tunneling in superfluid/superconducting systems is provided by the Feynman two-state 
model [69| : the macroscopic wavefunctions tp^ and ip R of the left and right systems obey the equations 



..dip L 



= E L ^ L - KtpR (Al) 

ih^=E R Tp R -Ki> L . (A2) 

The two-state model also describes also the tunneling of Bose-Einstein condensates in double well potentials [55]: the 
effect of the interactions between atoms in the wells results in cubic terms of the form [/|t/>l| ipL an d L/IV'itl ipR 
added to the right-hand sides of Eqs. ( |A1 |- dA2j ). In our case, since the D = 1 in it ial s tate p0| ) has mostly projections 
on the ground and first excited many-body states, we limit ourself to Eqs. (Al !-(A2| (with U = 0). 



Setting ip s 



(with s = L, R), the equations for z = (M L - M R )/jM L + Mr) and 



reads 



dz 
_dcf> 



-2Ky/l 
2Kz 
'dt ~ 7l == ? 



(A3) 
(A4) 



for the symmetric case El = E R . 
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The system (A3l-(A4| can be derived from the Hamiltonian 

H = -2Ky/l - z 2 cos0 
in which the time evolution of the conjugated variables <fr, z is found from 

dU 



hi 



h<t> 



dU 

dz 



By defining the angular variable 9 such that z = sin# G [—1,1], one finds from (A5) 



h6 = ^2Ksin<j) 
h<j) = 2K tan cos < 



(A5) 

(A6) 
(A7) 



(A8) 
(A9) 



where the =p sign accounts for the determination of the square root. The time-dependent relative occupation is a 
function of time only through the relative phase <f). Starting from ( A8 1-( A9 1 and identifying with a prime the derivative 
with respect to 0, one has 

M d6 d<j> n „ . l 



from which 



tan#- 



d9 



=p tan i 



By integration one obtains 



cos6> = 



A 



COS ( 



(A10) 



where the constant Aq — =p cos <fio cos #o is fixed by the initial conditions. Defining the current I as I — Ml one has 
/ = Mtz/2 where Mr = Ml + Mb is the total number of particles (pairs, in our case). Using (A3 1 one has 



TfA.\ MT h a KM T A 

Hp) = —^-ecose = — tan< 



(All) 



We conclude the Appendix by observing that for the linear two-state model here considered (U = 0) the phase 
difference does not overcome the value ir/2 (more precisely if \<fi(t — 0)| < 7r/2, then \4>(t)\ < 7r/2). 
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